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Abstract 

Probabilistic Boolean network (PBN) modelling is a semi-quantitative approach widely used for the study of the 
topology and dynamic aspects of biological systems. The combined use of rule-based representation and probability 
makes PBN appealing for large-scale modelling of biological networks where degrees of uncertainty need to be 
considered. 

A considerable expansion of our knowledge in the field of theoretical research on PBN can be observed over the past 
few years, with a focus on network inference, network intervention and control. With respect to areas of applications, 
PBN is mainly used for the study of gene regulatory networks though with an increasing emergence in signal 
transduction, metabolic, and also physiological networks. At the same time, a number of computational tools, 
facilitating the modelling and analysis of PBNs, are continuously developed. 

A concise yet comprehensive review of the state-of-the-art on PBN modelling is offered in this article, including a 
comparative discussion on PBN versus similar models with respect to concepts and biomedical applications. Due to 
their many advantages, we consider PBN to stand as a suitable modelling framework for the description and analysis 
of complex biological systems, ranging from molecular to physiological levels. 
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Background 

A large number of formal representation types that exist 
in Systems Biology are used to construct distinctive math- 
ematical models, each with their own strengths and 
weaknesses. On one hand, deciphering the complexity 
of biological systems by quantitative methods, such as 
ordinary differential equation (ODE) based mathemat- 
ical models, yields detailed representations with high 
predictive power. Such an approach is however often 
hampered by the low availability and/or identifiability 
of kinetic parameters and experimental data [1]. These 
limitations often result in the generation of relatively 
small quantitative network models. On the other hand, 
qualitative modelling frameworks such as the Boolean 
Networks (BNs), allow for describing large biological net- 
works while still preserving important properties of the 
systems [2]. The models pertaining to this latter class 
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fail nevertheless to offer a quantitative determination of 
the systems dynamics due to their inherent qualitative 
nature. 

Probabilistic Boolean networks (PBNs) were introduced 
in 2002 by Shmulevich et al. as an extension of the 
Boolean Network concept and as an alternative for mod- 
elling gene regulatory networks [3]. PBNs combine the 
rule-based modelling of a BN, as introduced by Kauff- 
man [4-7], with uncertainty principles, e.g., as described 
by a Markov chain [8]. In terms of applications, anal- 
ogously to the case of traditional BNs, the qualitative 
nature of state and time in a PBN framework allows 
for modelling of large-scale networks. The integrated 
stochastic properties of PBNs additionally enable semi- 
quantitative properties to be extracted. Existing analytic 
methods on PBNs allow for gaining a better under- 
standing of how biological systems behave, and offer 
in addition the means to compare to traditional BNs. 
Examples are the calculation of influences which rep- 
resent the quantitative strength of interaction between 
certain genes [3], or the determination of steady-state 
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distributions to quantitatively predict the activity of cer- 
tain genes in steady state [8] . 

It has been shown in the past years that the use of 
PBNs in the biological field is not limited to the molecu- 
lar level, but also can potentially be linked to applications 
in clinic. To name a few, Tay et al. constructed a PBN 
to demonstrate the interplay between dengue virus and 
different cytokines which mediate the course of disease 
in dengue haemorrhagic fever (DHF) [9]. Ma et al. pro- 
cessed functional Magnetic Resonance Imaging (fMRI) 
signals to infer a brain connectivity network comparing 
between Parkinsons disease patients and healthy subjects 
[10]. Even though the research efforts on PBNs in this 
direction are just sprouting, the results from such PBN 
studies can provide a first clue on a diseases etiology and 
progression. As PBNs are highly flexible for data integra- 
tion and as there exist a number of computational tools 
for PBN analysis, PBN is a suitable modelling approach 
to integrate information and derive knowledge from omic 
scale data which should in turn facilitate a physicians 
decision-making process in clinic. 

For the past decade, PBNs were the object of extensive 
studies, both theoretical and applied. Among theoretical 
topics, there are steady-state distribution, e.g., [11-13], 
network construction and inference, e.g., [14-16], net- 
work intervention and control, e.g., [17-19]. Several minor 
topics were investigated as well, including reachability 
analysis [20] or sensitivity analysis [21]. Other studies 
dealt with PBNs in biological systems at multi-level such 
as gene regulatory networks [22-24], signal transduction 
networks [25], metabolic networks [26], and also physi- 
ological networks [9,10] which could potentially link to 
medicine as previously mentioned. In parallel, a number 
of computational tools which facilitate the modelling and 
analysis of PBNs are also continuously developed [27-29]. 
Given the continuous development in this area due to 
the broad on-going range of research on PBNs, we offer 
a state-of-the-art overview on this modelling framework. 
A comparison of PBN to other graphical probabilistic 
modelling approaches is also enclosed, specifically with 
respect to Bayesian networks. Last but not least, a view 
of the theoretical and applied research on PBNs as mod- 
els for the study of multi-level biomedical networks is 
included. 

In order to provide a coherent overview of the recent 
advances on PBN, we start with several theoretical 
aspects, organised as follows: an introduction to PBNs and 
associated dynamics are given in Section Introduction to 
probabilistic Boolean networks and their dynamics', the 
construction and inference of PBNs as models for gene 
regulatory networks are presented in Section 'Construc- 
tion and inference of PBNs as models of gene regulatory 
networks', structural intervention and external control are 
discussed in Section 'Structural intervention and con- 



trol of PBNs', ending with the relationship between PBNs 
and other probabilistic graphical models in Section 'Rela- 
tionship between PBNs and other probabilistic graphical 
models! Later, in Section 'PBN applications in biological 
and biomedical studies' we present a broad summary of 
PBN applications as a representation of biological net- 
works followed by a discussion on the future applications 
of PBN in Systems Biology and Systems Biomedicine. A 
short conclusion is given in Section 'Conclusion'. 

Introduction to probabilistic Boolean networks 
and their dynamics 

Boolean networks 

A Boolean Network (BN) G(V, F), as originally introduced 
by Kauffman [4-7], is defined as a set of binary- valued 
variables (nodes) V = { X\f X^) . . . , OCyi I and a vector of 
Boolean functions/ = (fi,...,f n ). At each updating 
epoch, referred to as time point t (t = 0, 1, 2, . . .), the 
state of the network is defined by the vector x(t) = 
(x\ (t), X2(t), . . ., x n (t)), where Xi(t) is the value of variable 
Xi at time t, i.e., Xi(t) e {0, 1} (i = 1, 2, ... , n). For each 
variable X{ there exists a predictor set {Xi lf Xi 2 , . . ->Xi k(i) } 
and a Boolean predictor function (or simply predictor) f 
being the i-th element off that determines the value of X{ 
at the next time point, i.e., 

Xi(t + 1) =Mx h (t),x i2 (t), . . . ,*fc (0 (f)), (1) 

where 1 < i\ < i<i < • • • < i^ < n. Since 
the predictor functions of / are time-homogenous, the 
notation can be simplified by writing /(a?^, Xi 2 , . . . ,Xi k{i) )< 
Without loss of generality, k(i) can be defined to be 
a constant equal to n for all i by introducing ficti- 
tious variables in each function: the variable X{ is ficti- 
tious for a function / ]ff(xi, . . . , Xi-i, 0, . . . ,x n ) = 
f(x\, . . .,Xi-\, . . . ,x n ) for all possible values of 

xi, . . . , ...,x n .A variable that is not fictitious is 

referred to as essential The k(i) elements of the predictor 
set {Xi lf Xi 2 , . . . ,Xi k{f) } are referred to as the essential pre- 
dictors of variable Xi. The vector/ of predictor functions 
constitutes the network transition function (or simply the 
network function). The network function/ determines the 
time evolution of the states of the Boolean network, i.e., 
x(t + 1) = f(x(f)). Thus, the BN's dynamics is determin- 
istic. The only potential uncertainty is in the selection of 
the initial starting state of the network. 

Given an initial state, within a finite number of steps, 
the BN will transition into a fixed state or a set of states 
through which it will repeatedly cycle forever. In the first 
case, each such fixed state is called a singleton attractor, 
whereas in the second case, the set of states is referred to 
as a cyclic attractor. An attractor is either a singleton or 
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a cyclic attractor. The number of transitions required to 
return to a given state in an attractor is the cycle length of 
that attractor. The attractor structure of the BN is deter- 
mined by the particular combination of singleton and 
cyclic attractors, and by the cycle lengths of the cyclic 
attractors. The states within an attractor are called attrac- 
tor states, Non-attractor states are called transient and are 
visited at most once on any network trajectory. The states 
that lead into an attractor constitute its basin of attrac- 
tion. The basins form a partition of the state space of the 
BN. For example, in Figure 1 the state transition diagrams 
of four different Boolean networks with three variables 
are given (in fact all these Boolean networks constitute a 
probabilistic Boolean network — the framework of prob- 
abilistic Boolean networks is presented in Section '5'). For 
each of these networks attractor states and transient states 



are indicated and the cyclic- and singleton attractors are 
given. 

A Boolean Network with perturbations (BNp) is a BN 
with an introduced positive probability for which, at any 
transition, the network can depart from its current tra- 
jectory into a randomly chosen state, which becomes an 
initial state of a new trajectory. Formally, the perturba- 
tion mechanism is modelled by introducing a parameter 
p, 0 < p < 1, and a so-called perturbation vector y = 
(Kb Y2> • • • j Yn)> where yi> Yi*--->Yn are independent and 
identically distributed (i.i.d.) binary-valued random vari- 
ables a such that Pr{yi = 1} = p, and Pr{/i = 0} = 1 — p, 
for all i = 1, 2, . . . , n. For every transition step of the net- 
work a new realisation of the perturbation vector is given. 
If x(t) e {0, l} n is the state of the network at time t, then 
the next state x(t + 1) is given by either f(x(t)) or by 




(c) 



HOO.i 



:010' 



(d) 



MOOl 



'010i 



:oiii- 




1101 1 



■000r 



'01 1t- 



■>(boj) 



•101 J 



Figure 1 State transition diagrams of the four constituent Boolean networks of the PBN in Figure 2. For each constituent BN the attractor 
states and the transitions between them are indicated with solid circles and arrows, respectively. The remaining transitions and transient states are 
indicated with dashed arrows and circles, respectively, (a) The constituent BN of the PBN in Figure 2corresponding to transition function f ] . There is 
only one attractor, i.e., {01 1, 1 1 1}, which is a cyclic attractor. (b) The constituent BN of the PBN in Figure 2 corresponding to transition function f 2 . 
There are two cyclic attractors: {01 1, 1 1 1}, {001, 101} and one singleton attractor: {1 10}. (c) The constituent BN of the PBN in Figure 2 corresponding 
to transition function f 3 . {001 , 1 1 0, 1 1 1 } is the cyclic attractor. (d) The constituent BN of the PBN in Figure 2corresponding to transition function f 4 . 
There are two attractors: a cyclic one, i.e., {001 ,111} and a singleton one, i.e., {110}. 
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x(t) 0 y (t), where 0 is component-wise addition modulo 
2 and y(t) e {0, l} n is the realisation of the perturbation 
vector for the current transition. The choice of the state 
transition rule depends on the current realisation of the 
perturbation vector. Two cases are distinguished: either 
y(t) = 0 or at least one component of y(t) is 1, i.e., 
y (t) 7^ 0. In the first case, which happens with probability 
(1 — p) n , the next state is given by f(x(t)). In the second 
case, given with probability 1 — (1 — p) n , the next state 
is determined as x(t) 0 y(t): if yi = h then Xi changes 
its value; otherwise it does not (i = 1, 2, ... , n). Since 
y (t) 7^ 0, at least one of the nodes flips its value. 

The attractors of a Boolean network characterise its 
long-run behaviour [8] . However, if random perturbations 
are incorporated, the network can escape the attractors. 
In particular, perturbations allow the system to reach 
any of its states from any current state in one transi- 
tion. In consequence, the dynamics of the BNp is given 
by an ergodic Markov chain [30], b having a unique sta- 
tionary distribution which simultaneously is its steady- 
state (limiting) distribution. The steady-state probability 
distribution, where each state is assigned a non-zero 
probability, characterises the long-run behaviour of the 
BNp. Nevertheless, if perturbation probability is very 
small, the network will remain in the attractors of the orig- 
inal network for most of the time, meaning that attractor 
states will carry most of the steady-state probability mass 
[8]. In this way the attractor states remain significant for 
the description of the long-run behaviour of a Boolean 
network after adding perturbations. Thus, a BNp inherits 
the attractor-basin structure from the original BN; how- 
ever, once an attractor has been reached, the network 
remains in it until a perturbation occurs that throws the 
network out of it [31]. 

Probabilistic Boolean networks 

PBNs were introduced in order to overcome the deter- 
ministic rigidity of BNs [3,32,33], originally as a model for 
gene regulatory networks. A PBN consists of a finite col- 
lection of BNs, each defined by a fixed network function, 
and a probability distribution that governs the switching 
between these BNs. 

Formally, a probabilistic Boolean network G(V, F) is 
defined by a set of binary- valued variables (nodes) 0 V = 
{x\,X2, . • • ,x n } and a list of sets T = (F\,F2, • • • >Fn)< For 
i = 1,2,..., n the set F t is given as {/ 1 (0 ,/ 2 (0 , . . . 

where f^ l \ 1 < j < l(i), is a possible Boolean predictor 
function for the variable xu with /(/) the number of pos- 
sible predictors for X{. In general, each node X[ can have 
/(/) different sets of essential predictors, each specified for 
a particular predictor function in F;. A realisation of the 
PBN at a given instant of time is determined by a vec- 
tor of predictor functions, where the /th element of that 



vector contains the function selected at that time point 
for x^ For a PBN with N realisations there are N possible 
network transition functions fi,f 2 > • • • >/n °f tne f° rm 

// = (C4 2) '---'//f>' / = 1 > 2 '---' N ' i < h < m. 

e Fj, and ; = 1, 2, . . . , n. Each network function f l 
defines a constituent Boolean network, or context, of the 
PBN. 

Let/ = (f (1 \f {2 \ . . . ,f {n) ) be a random vector taking 
values in Fi x F 2 x • • • x F n ; in other words,/ is a random 
vector that acquires as value any of the realisations of the 
PBN. The probability that the predictor tf \ 1 < ; < /(/), 
is selected to determine the value of X{ is given by 

c »=Pr{/©=yf }= J2 K{f=fl)- (2) 

It follows that J2f=i cf = 1. The PBN is said to be 
independent if the random variables f^\f^ 2 \ . . . ,f^ are 
independent. Assuming independence, there are N = 
nf=i realisations (constituent BNs) of the PBN and the 
probability distribution on/ governing the selection of a 
particular realisation is given by Vv{f = f t } = nf=i - 
An example of a PBN with three nodes is given in 
Figure2. 

At each time point of the PBNs evolution, a decision 
is made whether to switch the constituent network. This 
is modelled with a binary random variable f: if § = 
0, then the current constituent network is preserved; if 
f = 1, then a context is randomly selected from all the 
constituent networks in accordance with the probability 
distribution of/. Notice that this definition implies that 
there are two mutually exclusive ways in which the context 
may remain unchanged: 1) either § = 0 or 2) § = 1 and 
the current network is reselected. The functional switch- 
ing probability q = Pr(§ = 1) is a system parameter. Two 
cases are distinguished in the literature: if q = 1, then 
a switch is made at each updating epoch; if q < 1, then 
the PBNs evolution in consecutive time points proceeds 
in accordance with a given constituent BN until the ran- 
dom variable § calls for a switch. If q = 1, as originally 
introduced in [32], the PBN is said to be instantaneously 
random; if q < 1, it is said to be context-sensitive. The 
former models uncertainty in model selection, the lat- 
ter models the situation where the model is affected by 
latent variables outside the model [34] . As an example let 
us consider the PBN given in Figure 2. Let the PBN be 
instantaneously random, i.e., q = 1. The four constituent 
BNs associated with the four transition functions /^/^ 
/ 3 , and/ 4 , are given in Figure 1. Further, let us assume 
that the initial state is the state 101 and that the con- 
secutive realisations ^ef v f 2 ,f 4 ,f 3 ,f 2 ,f 2 ,f 3 ,f 4 ,f 4 ,.... 
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Figure 2 An example of truth table, state transition diagram, and transition probability matrix of a PBN. The truth table, the state transition 
diagram, and the transition probability matrix A of a PBN without perturbations consisting of three variables V = {xi,X2,X3} and T = (F^.Fj.Fi), 



f(2) f( 2 )l 



Since there is one predictor function for node X] and two predictors for nodes x 2 and x$, 



(f 0) /f (2) /f (3) )//2 = (f; 



where F ] = {f| lj },F 2 = {f\ l) ' f 2 

there are 1 -2-2 = 4 realisations of the PBN given by four network transition functions f ] 

f 3 = (fl ]) J^J^), and f 4 = (f^\ f 2 (2) , f 2 (3) ) with associated probabilities q =0.1 2, c 2 =0.1 8, c 3 =0.28, and c 4 = 042, respectively. For example, 
= 1 • 0.7 • 0.4 = 0.28. The edges in the state transition diagram are labelled with the transition probabilities. As can be seen from 



(1) f (2) 



(1) J2) (3) 



C 3 =C ] -l 2 - L1 

the state transition diagram, the underlying Markov chain is irreducible and aperiodic, thus ergodic. The steady-state (limiting) distribution for the 
chosen q values, /' = 1 A, is given by [ j^g, f^fy, -^j, ^§7, f^W\^ ^ e states are considered in the lexicographical order from 000 

to 111). ' 



Then, the corresponding time evolution of the PBN (tra- 
jectory) is given by the following sequence of state tran- 
sitions: 101 -> 001 -> 110 -> 110 -> 111 -> 011 -> 
111 -> 001 -> 100 -> 011 -> .... Irrespective of which 
constituent network (realisation) is selected next, the con- 
secutive state in the trajectory is going to be 111 as the 
probability of moving from 011 to 111 is C1+C2+C3+C4 = 1. 

A Probabilistic Boolean Network with perturbations 
(PBNp) is the variant of the PBN framework in which 
each constituent network is a BNp with a common per- 
turbation probability parameter p, 0 < p < 1, and a 
perturbation vector y. If x(t) e {0, l} n is the current state 
of the network and y (t) = 0, then the next state of the 
network is determined according to the current network 
function f [} i.e., %{t + 1) = //(*(*)). If x(t) e {0, \} n is 
the current state and y(t) ^ 0, then x(t + 1) = x(t) 0 
y(t). Whereas a context switch in a PBNp corresponds 
to a change in latent variables, resulting in a structural 
change in the functions that govern the PBNp, a random 
perturbation reflects a transient value change that leaves 
the network wiring unmodified, as for example in the 
case of gene activation or inactivation caused by external 
stimuli such as stress conditions or small molecule 
inhibitors [8]. 



The relationship between the four frameworks, i.e., 
Boolean networks, Boolean networks with perturbations, 
probabilistic Boolean networks, and probabilistic Boolean 
networks with perturbations is schematically depicted in 
Figure 3. 

Dynamics of PBNs 

A Boolean network with perturbations can be viewed as 
a homogenous irreducible Markov chain X t , with state 
space X = {0, 1} W , where n is the number of nodes in the 
BNp. Let P y (x) = Pr[X; 0+ i = x\X to = y] be the Markov 
chain transition probability from state y to state x at any 
instant to. This probability is a weighted sum of two tran- 
sition probabilities, one for the BN, with probability (1 — 
p) n , and the other for the perturbations, with probability 
l-(l-#,i.e, 

P,(*) = l[f<y)^](l-/0*Kl-l|^ 

(3) 

where p is the perturbation probability, 1 is the indicator 
function (lrp] = 1 if the proposition P is true, and l[p] = 0 
otherwise), and rj(x,y) is the Hamming distance between 
the binary vectors x and y. 
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Figure 3 Relationships between the frameworks of Boolean and 
probabilistic Boolean networks. A Boolean network (BN) can be 
converted to a Boolean network with perturbations (BNp) by 
introducing a probability parameter p, 0 < p < 1 , and a perturbation 
vector (y). A probabilistic Boolean network (PBN) is built upon a 
number of constituent BNs and a probability distribution governing 
the choice of the Boolean network in accordance with which the next 
transition is made. Analogically, a PBN can be converted to a 
probabilistic Boolean network with perturbations (PBNp) by 
introducing a probability parameter p, 0 < p < 1 , and a perturbation 
vector (y). A probabilistic Boolean network (PBN) is built upon a 
number of constituent BNps and a probability distribution governing 
the choice of the BNp in accordance with which the next transition is 
made. 



The Markov chain X t is ergodic, which follows from 
the fact that it is aperiodic, irreducible, and defined on 
a finite state space. In other words, it possesses a unique 
stationary distribution, being simultaneously its steady- 
state (limiting) distribution. If Vy\x) is the probability 
that the system transitions from y to x in t time steps, i.e., 
Vy\x) = Vv[X tQ+t = x\X to = y], then the steady-state 
distribution n oiX t is defined by 7t(x) = lim^oo Pj?(#) 
for any initial state k e X. For a set of states B c X the 
steady-state probability is given by tt(B) = J^xeB 71 ^ ~ 
lim^ooPf (B) for any initial state k e X. For exam- 
ple, the steady- state distribution of the Markov chain 
given by the transition probability matrix in Figure 2 is 

r _7_ 3640 49 716 175 238 2548 4696 i (fh ctatpc 
L 1609' 14481' 4827' 4827' 4827' 4827' 14481' 14481 J vulc * LcLLC * 

are considered in the lexicographical order from 000 to 
111). 

In the case of a probabilistic Boolean network, the tran- 
sition probabilities P y (x) of the underlying Markov chain 
Xt depend on the probability of selecting a network tran- 
sition function k = 1, 2, . . .,N, that determines the 
transition from y to x i.e., 

N 

V y (x) = Pr[X m = x\X t =y] = J2 l \f k (y)=*Y^\f =fh)> 

k=l 

(4) 



where N, as before, is the number of constituent BNs and/ 
is a random vector determining the PBNs realisation. Let- 
ting x and y range all states in X } the transition probability 
matrix A of size 2 n x 2 n can be formed and expressed as 

N 

A = J2^k'V*{f=f k }> (5) 

k=l 

where is the transition matrix corresponding to the 
/c-th constituent BN. 

Now, adding perturbations with probability p makes 
the underlying finite-space Markov chain X t of the PBNp 
aperiodic and irreducible, hence ergodic. This allows the 
network dynamics of a PBNp to be studied with the use 
of the rich theory of ergodic Markov chains [30]. In par- 
ticular, in the case of instantaneously random PBNps, the 
transition probability matrix A is given by 

A= (l-pf-A + P, (6) 

where P is the perturbation matrix of the form 

p y>x = (i - i [x=y] ) P ^y\i - p f-^y), ( 7 ) 

where, as before, 1 is the indicator function and r) is the 
Hamming distance. As in the case of BNps, the ergod- 
icy of the underlying Markov chain ensures the existence 
of the unique stationary distribution being the limiting 
distribution of the chain. 

By definition, the set of attractors of a PBN is the union 
of the sets of attractors of the constituent networks [8]. 
Notice that whereas in a BN two attractors cannot inter- 
sect, attractors from different contexts can intersect in 
the case of a PBN. Similarly as in the case of Boolean 
networks, attractors play a major role in the characterisa- 
tion of the long-run behaviour of a probabilistic Boolean 
network. If, however, perturbations are incorporated, the 
long-run behaviour of the network is characterised by its 
steady-state distribution. Nevertheless, if both the switch- 
ing and perturbation probabilities are very small, then the 
attractors still carry most of the steady-state probability 
mass [8]. From a biological point of view attractors of such 
networks are interesting as they can be given a clear bio- 
logical interpretation: they can be used to model cellular 
states [31]. For example, in the context of gene regulatory 
networks, it is believed that attractors can be interpreted 
as cellular phenotypes [7,8]. Thus, the long-run behaviour 
of the network given by its steady-state probabilities is 
of a special interest. Specifically, the attractor steady- 
state probabilities, i.e., 7t(A), where A is an attractor, are 
important. There are a number of approaches towards the 
determination and analysis of the steady-state distribution 
of a PBNp. We review them shortly. 

First, one approach to the steady-state analysis is to con- 
struct the state transition matrix in some form or another 
and then apply some numerical methods, e.g., iterative, 
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decompositional or projection methods [35]. A transi- 
tion matrix based approach in which the sparse transition 
matrix is constructed in an efficient way and the so- 
called power method, which is applied to compute the 
steady-state probability distribution, is proposed in [36]. 
Unfortunately, the size of the state space grows expo- 
nentially in the number of nodes (genes) and becomes 
prohibitive for matrix-based numerical analysis of larger 
networks [11]. In [12], an approximation method for com- 
puting the steady-state probability distribution of a PBNp 
is derived from the approach of [36]. This method neglects 
some constituent BNps with very small probabilities dur- 
ing the construction of the transition probability matrix. 
An error analysis is given to demonstrate the effective- 
ness of this approach. Further, in [13] and [37] a matrix 
perturbation method for computing the steady-state 
probability distribution of PBNps is proposed together 
with its approximation variant. The proposed meth- 
ods make use of certain properties of the perturbation 
matrix, P. 

Second, Markov chain Monte Carlo methods [38] rep- 
resent a feasible alternative to numerical matrix-based 
methods for obtaining steady-state distributions. Given an 
ergodic Markov chain, a Monte Carlo simulation method 
has been proposed: the probability of being in state x in 
the long run can be estimated empirically by simulating 
the network for a sufficiently long time and by count- 
ing the percentage of time the chain spends in that state 
regardless of the starting state [8]. A set of examples of 
Monte Carlo simulations from the PBN example in Figure 
2 is shown in Figure 4. However, the question that remains 
is how to judge whether the simulation time is sufficiently 
long? The key factor here is the convergence, which in the 
case of a PBNp is known to depend to a large extent on 
the perturbation probability p [11]. Several approaches for 
determining the number of iterations necessary to achieve 
convergence were developed. A typical class consists of 
methods based on the second-largest eigenvalue of the 
transitions probability matrix, but due to reasons already 
mentioned above, these approaches can be impractical 
for larger networks. Another method utilises the so-called 
minorisation condition for Markov chains [39] to provide 
a priori bounds on the number of iterations. However, the 
usefulness of this approach is also limited (see [11] for 
details). There exist a number of methods for empirically 
diagnosing convergence to the steady-state distribution 
[40,41]. In [11] two of them are considered: one, based 
on the Kolmogorov-Smirnov test, a nonparametric test 
for the equality of continuous, one-dimensional proba- 
bility distributions, and, second, the approach proposed 
in [42] which reduces the study of convergence of the 
chain to the investigation of the convergence of a two- 
state Markov chain. For illustration of application of these 
approaches to PBNs, we refer to [11] where the joint 



steady-state probabilities of combinations between two 
genes in human glioma gene expression data set were 
analysed. 

Finally, as shown in [31], analytical expressions for the 
attractor steady-state probabilities can be derived both 
for BNps and PBNps. The obtained formulas are fur- 
ther exploited to propose an approximate steady-state 
computation algorithm. 

We just shortly mention here that in the case of 
probabilistic Boolean networks without perturbations the 
dynamics is given by a Markov chain that does not nec- 
essarily be ergodic, specifically the Markov chain may 
contain more than one so-called ergodic set of states, also 
referred to as a closed, irreducible set of states in the lit- 
erature. An ergodic set of states C in a Markov chain 
is defined as a set of states where all states communi- 
cate and no state outside C is reachable from any state 
in C d . The notion of an ergodic set of the correspond- 
ing Markov chain in probabilistic Boolean networks is the 
stochastic analogue of the notion of an attractor in stan- 
dard Boolean networks [32]. Notice, however, that the 
ergodic sets and the attractors of a PBN or PBNp may dif- 
fer. In the case of probabilistic Boolean networks without 
perturbations where the underlying Markov chain con- 
tains more than one ergodic set, considering the ergodic 
sets rather than the attractors may be more significant 
for understanding the long-run behaviour of the net- 
work. For example, in the context of modelling biolog- 
ical processes with PBNs, cellular phenotypes may in 
fact be represented by the ergodic sets. For more details 
see [32,43,44]. 

A number of other issues related to probabilistic 
Boolean network dynamics have been considered in the 
literature. We briefly list them here. In [45,46], the 
ordering of network switching and state transitions in 
context-sensitive PBNs are considered and its influence on 
the steady-state probability distributions is investigated. 
Algorithms for enumeration of attractors in probabilistic 
Boolean networks are discussed in [47]. Stability and sta- 
bilisation issues of PBNs are covered in [48] . Further, net- 
work transformations from one to another without losing 
some crucial properties, e.g., the steady-state probability 
distribution, are considered in [49]. For this purpose the 
concepts of homomorphisms and 6-homomorphisms for 
probabilistic regulatory networks, in particular PBNs, are 
developed. 

Construction and inference of PBNs as models of 
gene regulatory networks 

One approach to the dynamical modelling of gene regula- 
tion is based on the construction and analysis of network 
models. Generally, in the study of dynamical systems, 
long-run behaviour characteristics are of utter impor- 
tance and their determination is a main aspect of system 



Trairatphisan et al. Cell Communication and Signaling 201 3, 1 1 :46 
http://www.biosignaling.eom/content/1 1/1/46 



Page 8 of 25 



(a) 



1 

0.8 
0.6 
0.4 
0.2 
0 



(C) 



0.8 
c\j 0.6 

X 

0.4 
0.2 



mean at steady state = 1 



5 

time 



10 



mean at steady state = 0.66 



5 

time 



10 



(b) 



1 

0.8 
0.6 
0.4 
0.2 
0 



mean at steady state = 0.5 



5 
time 



10 




Figure 4 Dynamical simulations of node xi of the example network in Figure 2, with initial state k = 000. (a) Dynamics of xj governed by 
the constituent BN corresponding to the transition function f ] , where C] = 1 , Cj = C3 = C4 = 0. Starting from 000 the periodic attractor {011,111} 
is reached. The probability of {xj = 1 } given by the stationary distribution is 1 . (b) Dynamics of xj governed by the constituent BN corresponding to 
the transition function f 4 , where c 4 = 1 , Q = c 2 = c 3 = 0. Starting from 000 the periodic attractor {001 , 1 1 1 } is reached. The probability of {x 2 = 1 } 
given by the stationary distribution related to the reached attractor, i.e., [0, ^,0, 0, 0,0, 0, \] (the states are considered in the lexicographical order), is 
0.5. (c,d) Examples of X2 dynamics in the full PBN as given in Figure 2. Starting from 000 different trajectories are obtained for different simulation 
runs. The underlying Markov chain is ergodic and a unique stationary distribution, being the steady state (limiting) distribution, exists therefore. The 
steady state probability of {xj = 1 } is 0.66. 



analysis. Reversely, the task of constructing a network 
possessing a specific set of properties is a subject of sys- 
tem synthesis. However, this inverse problem is usually 
ill-posed, i.e., there may be many models, or none, with 
the given properties [50]. Here we concentrate on the 
problem of inference from data in the framework of prob- 
abilistic Boolean networks, an inverse problem in which 
a network is constructed relative to some relationship 
with the available data. An outline of the workflow in 
network inference in the PBN framework is shown in 
Figure 5. 

A data-driven approach for model construction con- 
sists of inferring the model structure and model param- 
eters from measurement data, which in the case of gene 
regulation most commonly are gene expression measure- 
ments obtained with microarray technology. However, 
such data are continuous in nature. Thus, prior to the 
inference of Boolean or other discrete-type models (e.g., 
ternary) the measurements are usually discretised. The 
most common discretisation is binary (0 or 1) or ternary 
(usually -1, 0, 1) [8]. Discretisation is often justified as 
biological systems commonly exhibit switch-like on/off 



behaviour. Moreover, there are also a number of prag- 
matic reasons for quantising the measurements, e.g., it 
reduces the level of model complexity implying less com- 
putation and lower data requirements for model identi- 
fication, provides a certain level of robustness to noise 
in the data, and has been shown to substantially reduce 
error rates in microarray-based classification [8,51-53]. A 
number of methods for discretisation of gene expression 
data exist, many of them having their origin in signal pro- 
cessing. One approach to quantisation was proposed in 
[54]: given some thresholds X\ < T2 < ... (e.g., cor- 
responding to limiting cases of a sigmoidal response), a 
multilevel discrete variable x is defined as x = cp(x) = r^ 
for < x < T£ + i. As mentioned in [8], the thresh- 
olds can either come from prior knowledge or be chosen 
automatically from the data. In fact, there are various 
ways for optimal selection of the thresholds t>. One of 
the most popular methods is the Lloyd-Max quantizer, 
which amounts to minimising a so-called mean square 
quantisation error, see [55] for details. Approaches spe- 
cific to binarising gene expression data can be found 
in [56-58]. Recently, Hopfensitz et al. [58] proposed a 
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Figure 5 An outline of the workflow in network inference and control in the PBN framework. Microarray data, either from steady-state or 
time-course measurements, are typically binarised or discretised into discrete values. A heuristic approach, such as using genetic algorithms, is 
generally applied to identify constituent Boolean networks of the inferred PBN. Regularisation methods can be further applied to improve the 
accuracy of the inference with use of prior information on the network structure or dynamical rules. A number of well-established methods are 
subsequently applied to determine the predictor probability of each constituent Boolean network, thus the PBN is inferred. The inferred PBN can 
subsequently be perturbed with the methods on structural intervention or external control. The goal of network control is to increase the 
probability of reaching desirable states in the corresponding PBN. 



new approach to binarisation which incorporates mea- 
surements at multiple resolutions. The method, called 
Binarization across Multiple Scales, is based on the com- 
putation of a series of step functions, detection of the 
strongest discontinuity in each step function and the esti- 
mation of the location and variation of the strongest 
discontinuities. Two variants of the method are proposed 
which differ in the approach towards the calculation of 
the series of step functions. The proposed method allows 
thresholds determination even with limited number of 
samples and simultaneously provides a measure of thresh- 
old validity - the latter can further be used to restrict 



network inference only to measurements yielding rele- 
vant thresholds. An example of application of binarisation 
to real data in the context of modelling with PBNs can 
be found in [10], where a brain connectivity network of 
Parkinsons disease is analysed. Binarisation is performed 
on fMRI real-valued data along the method recently 
proposed in [59]. 

One of the most straightforward inferential approaches 
is the consistency problem (also referred to as the extension 
problem), that entails a search for a rule from experimental 
data [8,60-62]. The problem amounts to finding in a spec- 
ified class of Boolean functions one that complies with 
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two given sets of "true" and "false" Boolean vectors, i.e., 
a function that takes the value 1 for each of the "true" 
vectors and 0 for each of the "false" vectors. 

In the case of real experimental data, a consistent exten- 
sion may not exist either due to measurement noise or 
due to some underlying latent factors or other external 
influences not considered in the model [8]. In such case 
instead of searching for a consistent extension a Boolean 
function that minimises the number of misclassifications 
(errors) is considered. This problem is known as the 
best-fit extension problem [61] and is computationally 
more difficult than the consistency problem, since the 
latter is a special case of the former. 

The application of PBN for modelling of large-scale 
networks is often impeded by limited sample sizes of 
experimental data. As mentioned in [63], main challenges 
in automated network reconstruction arise from the expo- 
nential growth of possible model topologies for increasing 
network size, the high level of variability in measured 
data often characterised by low signal to noise ratios, and 
the usually large number of different components that 
are measured versus relatively small number of differ- 
ent observations under changing conditions, e.g., number 
of time points or perturbations of the biological system. 
Together these problems lead to non-identifiability and 
over- fitting of models [63]. In such cases any prior infor- 
mation on the network structure or dynamical rules is 
likely to improve the accuracy of the inference [8,64]. 
This information usually pertains to model complexity 
and is used to penalise excessively complex models. For 
this purpose, the so-called regularisation methods can 
be employed. The most popular regularisation assump- 
tion in gene regulatory modelling is that the inferred 
models should be sparse, i.e., the number of regulators 
acting on a gene is low [65-68] or that the node degree 
in biological networks is often power law distributed, 
with only few highly-connected genes, and most genes 
having small number of interaction partners [63,69]. Reg- 
ularisation is a well-established inference approach in the 
framework of Bayesian networks (see, e.g., [63,70,71]) and 
can be also used in the framework of BNs and PBNs. 
For example, in the case of inference of Boolean net- 
works, the so-called sensitivity regularisation method has 
been proposed [64]. Due to limited sets of data, the 
estimates of the errors of a given model in the best- 
fit extension problem, which themselves depend on the 
measurements, may be highly variable [64]. The regu- 
larisation is built on the observation that the expecta- 
tion of the state transition error generally depends on 
a number of terms, among others the sensitivity devi- 
ation which is a difference in the sensitivities of the 
original and the inferred networks. In consequence, as 
argued in [64], the sensitivity deviation can be incorpo- 
rated as an additional penalty term to the best- fit objective 



function, reflecting the hypothesis that the best inference 
should have a small error in both state transition and 
sensitivity. 

In order to infer a PBN, strong candidates for regu- 
lar Boolean networks need to be identified first. This 
can be performed with generic methods mentioned in 
[72] such as literature data compilation, the gene associ- 
ation networks approach [73,74] or by applying a heuris- 
tic approach, e.g., a genetic algorithm, which searches 
through the model space to find good candidates for 
the network structure with respect to a specified fitness 
function. Next, the candidates' predictor functions are 
combined into a set of network transition functions for the 
PBN. An example of PBN model selection using heuristics 
can be found in [75]. 

A common strategy for determining the predictor prob- 
abilities relies on the coefficient of determination (CoD) 
between target and predictor genes [8,32,72,76]. The CoD 
is a measure of relative decrease in error from estimat- 
ing transcriptional levels of a target gene via the levels 
of its predictor genes rather than the best possible pre- 
diction in the absence of predictor genes [8]. The CoDs 
can be then translated to the predictor probabilities. How- 
ever, as pointed out in [77], for each gene, the maximum 
number of possible predictors as well as the number of 
their corresponding probabilities is equal to 2 2 , where 
n is the number of nodes. This implies that the number 
of parameters in the PBN model is 0(n2 2 ) e . Therefore, 
the applicability of the CoD approach is significantly lim- 
ited due to the model complexity or imprecisions owing 
to insufficient data sample size. This hindrance is often 
surpassed by imposing some constraints on the maximum 
size of admissible predictors for each gene. 

In [50] the authors consider the attractor inverse prob- 
lem, that involves designing Boolean networks given 
attractor and connectivity information. Two algorithms 
for solving this problem are proposed. They are based 
on two assumptions on the biological reality: first, the 
biological stability, i.e., that most of the steady-state prob- 
ability mass is concentrated in the attractors and, second, 
the biological tendency to stably occupy a given state, 
i.e., attractors are singleton attractor cycles consisting of 
a single state. The first algorithm operates directly on the 
truth table, while taking into account simultaneously the 
information on the attractors and predictor sets. There is 
however no control on the level-set structure. The sec- 
ond algorithm works on the state transition diagram that 
satisfies the design requirements on attractor and level- 
set structures and checks whether the associated truth 
table has predictor sets that agree with the design goals. 
The proposed algorithms can be further used in a pro- 
cedure for designing PBN from data. In the approach 
described in [50], a collection of BNs is generated by 
the first algorithm, then some of the BNs are selected 
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based on the basin sizes criterion and combined in a 
PBN whose steady-state distribution closely matches the 
observed data frequency distribution. This design pro- 
cedure has been applied to gene-expression profiles in a 
study of 31 malignant melanoma samples in [50]. 

An inverse PBN construction approach is also described 
in [78]. This work relies on expressing the probabil- 
ity transition matrix as a weighted sum of Boolean 
network matrices. A heuristic algorithm with 0(m2 n ) 
complexity is proposed, where n, m stand for the 
number of genes, respectively the number of non-zero 
entries in the transition matrix. The authors also intro- 
duce an entropy based probabilistic extension, both 
algorithms being analysed against random transition 
matrices. 

Usually, the optimal predictor for a gene will not be 
perfect as there will be inconsistencies in the data. In 
[79] it is proposed to model these inconsistencies in 
a way that mimics context changes in genomic regu- 
lation, with the intention to view data inconsistencies 
as caused by latent variables. The inference procedure 
of [79] results in PBNs whose contexts model the data 
in such a way that they are consistent within each 
context. The key criterion for network design is that 
the distribution of data states agrees with the distribu- 
tion of expected long-term state observations for the 
system. 

The probabilities of the system being in a particular 
context and the number of constituent networks are deter- 
mined by the data. The approach of [79] can be seen as 
imposing a structure on a probabilistic Boolean network 
that resolves inconsistencies in the data arising from mix- 
ing of data from several contexts. It should be noted that 
in this approach the contexts are determined directly by 
the data, whereas in [32] and [80] constituent networks 
depend on the number of high-CoD predictor sets or 
high Bayes-score predictor sets, respectively, and these 
in turn depend on the designers choice of a threshold. 
Moreover, the number of constituent networks is deter- 
mined by how inconsistencies appear in the data, not 
the number of states appearing in the data (see [8] for 
an example). The contextual-design method of [79] has 
been applied to expression profiles for melanoma genetic 
network. 

We just mention here that also information theoretic 
approaches were considered for inference of PBN from 
data. Probably the most widely studied methods are based 
on the minimum description length (MDL) principle [81]. 
Descriptions of inference algorithms that utilise this prin- 
ciple can be found, e.g., in [8,82,83]. 

The manner of inference depends on the kind of exper- 
imental data available. There are two cases: 1) time-series 
data and 2) steady-state data. We proceed with presenting 
them briefly. 



Time-course measurements 

It is assumed that the available data are a single temporal 
sequence of network states. In this case, given a suffi- 
ciently long sequence of observations, the goal is to infer 
a PBN that is one of plausible candidates to have gener- 
ated the data. Usually, an inference procedure for this type 
of problem constructs a network that is to some extent 
consistent with the observed sequence. 

In [84,85], the inference in case of context-sensitive 
PBNs with perturbations is considered, where the proba- 
bility of switching from the current constituent Boolean 
network to a different one is assumed to be small. The 
proposed inference procedure consists of three main 
steps: first, identification of subsequences in the tempo- 
ral data sequence that correspond to constituent Boolean 
networks with use of so-called purity functions'; sec- 
ond, determination of essential predictors for each subse- 
quence by applying an inference procedure based on the 
transition counting matrix and a proposed cost function; 
finally, inference of perturbation, switching, and selec- 
tion probabilities. However, the amount of temporal data 
needed for inference with this approach is huge, especially 
due to the perturbation and switching probabilities: if they 
are very small, then long periods of time are needed to 
escape attractors and if they are large, estimation accu- 
racy is harmed. As stated in [85], if one does not wish to 
infer the perturbation, switching, and selection probabili- 
ties, then constituent-network connectivity can be discov- 
ered with decent accuracy for relatively small time-course 
sequences. 

A more practical way of inferring PBN parameters 
from time-course measurements is presented in [77]. The 
authors propose a multivariate Markov chain model to 
infer the genetic network, develop techniques for esti- 
mating the model parameters and provide an efficient 
method of estimating PBN parameters from their multi- 
variate Markov chain model. The proposed technique has 
been tested with synthetic data as well as applied to gene 
expression data of yeast. 

Further, in [86] the problem of PBN context estimation 
from time-course data is considered. The inference is con- 
sidered with respect to minimising both the conditional 
and unconditional mean-square error (MSE). The author 
proposes a novel state-space signal model for discrete- 
time Boolean dynamical systems, which includes as spe- 
cial cases distinct Boolean models, one of them being 
the PBN model. A Boolean Kalman Filter algorithm is 
employed to provide the optimal PBN context switch- 
ing inference procedure in accordance to minimisation of 
MSE. 

Steady-state data 

Here we consider a long-run inverse problem in the 
context of probabilistic Boolean networks as models for 



Trairatphisan et al. Cell Communication and Signaling 201 3, 1 1 :46 
http://www.biosignaling.eom/content/1 1 /I /46 



Page 12 of 25 



gene regulation. On one hand, in the case of microarray- 
based gene-expression studies it is often assumed that 
the data are obtained by sampling from a steady state. 
On the other hand, attractors represent the essential 
long-run behaviour of the modelled system [31]. Thus, 
in the modelling framework of Boolean networks it is 
expected that the observed data states are mostly the 
attractor states of a model network. In consequence, much 
of the steady-state distribution mass of the model net- 
work should lie in the states observed in the sample 
data [50,80,87]. In the case of Boolean networks with 
perturbations or probabilistic Boolean networks with per- 
turbations, the underlying dynamical system is an ergodic 
Markov chain, hence possesses a steady-state distribution. 
However, by imposing some mild stability constraints that 
reflect biological state stability, also in these frameworks 
most of the steady-state probability mass is carried by the 
attractors [31]. 

There are however inherent limitations to the con- 
struction of dynamical systems from steady-state data. 
Although the steady-state behaviour restricts the net- 
work dynamics, it does not determine the steady-state 
behaviour: there may be a collection of compatible net- 
works with a given attractor structure. In particular, it 
does not determine the Boolean networks basin structure. 
As a consequence, obtaining good inference relative to the 
attractor structure does not necessary entail valid infer- 
ence with respect to the steady-state distribution as the 
steady-state probabilities of attractor states depend on the 
basin structure [50,80]. In fact building a dynamical model 
from steady-state data is a kind of over- fitting [88]. 

Although the CoD has been used for inference of PBNs 
from steady-state data in [32], a fundamental problem is 
that the CoD cannot provide information on the direc- 
tion of prediction without time-course data. The resulting 
bidirectional relationships can affect the inferred graph 
topology by introducing spurious connections. Moreover, 
they can lead to inference of spurious attractor cycles that 
do not correspond to any biological state [8]. As a conse- 
quence, this suppressed the use of the CoD as a inference 
method for steady-state data. 

The inference methods that replaced the CoD approach 
are primarily based on the attractor structure [50,79] or 
graph topology [89]. In the former case, the key concern 
is to infer an attractor structure close to that of the true 
network. In the latter case, the focus is on the agree- 
ment between graph connections, e.g., as measured by 
the Hamming distance between the regulatory graphs [8] . 
In [16], an approach that achieves both preservation of 
attractor structure and connectivity based on strong gene 
prediction has been proposed. 

Another approach to the problem of constructing gene 
regulatory networks from expression data using the PBNs 
framework is proposed in [90]. The key element of this 



method is a non-linear regression technique based on 
reversible-jump Markov chain Monte Carlo (MCMC) 
annealing for predictor design. The network construc- 
tion algorithm consists of the following stages. First, for 
each target gene x\ (i = 1, 2, ... , n) in the network of 
n genes a collection of predictor sets is determined by 
applying a clustering technique based on mutual informa- 
tion minimisation. Optimisation f is performed with use 
of the simulated annealing procedure. This step reduces 
the class of different predictor functions available for 
each target gene. Next, each predictor set is used to 
model a predictor function fj^ by a perceptron con- 
sisting of both a linear and a nonlinear term, where 
k = 1, 2, ... , /(/), with /(/) the number of predictor sets 
found in the previous step for target gene %{. A reversible 
MCMC technique is used to calculate the model order 
and the parameters. Finally, the CoD is used to compute 
the probability of selecting different predictors for each 
gene. For a detailed description of this algorithm and its 
application to data on transcription levels in the context 
of investigating responsiveness to genotoxic stresses see 
[90]. It should be noticed that the proposed reversible- 
jump MCMC model for predictor design extends the 
binary nature of PBNs allowing for a more general model 
containing non-Boolean predictor functions that operate 
on variables with any finite number of possible discrete 
values [72]. 

As an alternative to the technique of [90], a fully 
Bayesian approach (without the use of CoD) for con- 
structing probabilistic gene regulatory networks, with an 
emphasis for network topology, is proposed in [80]. In 
this approach, the predictor sets of each target gene are 
computed, the corresponding predictors are determined, 
and the associated probabilities, based on the nonlinear 
perceptron model of [90], are calculated by relying on 
a reversible jump MCMC. Then, a MCMC method is 
used to search for the network configurations that max- 
imise the Bayesian scores to construct the network. As 
stated in [8], this method produces models whose steady- 
state distribution contains attractors that are either iden- 
tical or very similar to the states observed in the data. 
Moreover, many of the attractors are singleton attractors, 
which reflect the biological propensity to stably occupy 
a given state. The approach of [90] has been applied to 
gene-expression profiles resulting from the study of 31 
malignant melanoma samples presented in [91]. 

In [92] the inverse problem of constructing instanta- 
neously random PBNs from a given stationary distribution 
and a set of given Boolean networks is considered. Due 
to large size of this problem, it is formulated in terms of 
constrained least squares and a heuristic method based on 
Conjugate Gradient is proposed as a solution. 

In [93], the inverse problem of PBNs with perturba- 
tions is considered, where a modified Newton method is 
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proposed for computing the perturbation probability p 
where the transition probability matrix A and the steady- 
state probability of the PBNp x are known. The new 
algorithm makes use of certain properties of the set of 
steady-state nonlinear equations, i.e., Ax — x = 0, with 
p as the unknown variable. Considering these proper- 
ties improves the computational efficiency with respect 
to a direct approach in which every of the 2 n equations 
(n being the number of nodes) is solved and common 
solutions are reported. 

Structural intervention and control of PBNs 

Using PBNs for the modelling and analysis of biological 
systems can lead to a deeper understanding of the dynam- 
ics and behaviour of these systems (see Section 'Dynamics 
of PBNs'), paving the way for different methods used for 
system structure inference and data measurement (see 
Section 'Construction and inference of PBNs as models 
of gene regulatory networks'). Another major objective of 
such studies is to predict the effect a perturbation or an 
intervention has on the system structure, e.g., allowing to 
identify potential targets for therapeutic intervention in 
diseases such as cancer. Intervention strategies in PBNs, 
e.g., as to change the long-run behavior of networks in 
order to decrease the probability of entering some unde- 
sired state, rely on two different kinds of direction - 
structural intervention [8,33] and external control [8,18]. 
While the first approach can alter the underlying network 
structure permanently, the second one uses external con- 
trol to modulate the network dynamics. A classification of 
network control methods in the PBN framework is shown 
in Figure 5. 

Structural intervention 

The problem of performing a structural intervention in 
a PBN looks at how the steady-state probability of cer- 
tain states can be changed with only minimal structural 
modifications [8,33]. A more formal description is offered 
in the following. Given a PBN and two subsets A and 
B of its states, the associated steady-state probabilities 
tt(A), 71(B), have to be modified such as to approach some 
given values Xa, respectively A.#. This can be achieved by 
replacing the predictor function^ (of gene i in context 
k) with a new function gu a while keeping all other net- 
work parameters unchanged. We denote the steady-state 
distribution of the resulting PBN as /x. Then, it is possi- 
ble to interpret the problem as an optimisation one: given 
the state sets A, B, and two values Xa > 0, Xb > 0, 
with Xa + Xb < 1, find a context /<, a gene /, and a func- 
tion gfc to replace^, such as to minimises e(A,B) =\ 
/jl(A) — Xa I + | /jl(B) — Xb I, with respect to all contexts, 
genes, and predictor functions. Note that A and B can 
be used to represent both desirable as well as undesirable 
states. While this approach allows changing one predictor 



function at a time, a generalisation can be made by allow- 
ing a number of predictor functions or by adding more 
constraints on the selected functions, only to give a few 
examples. 

Shmulevich et al. [33] proposed using genetic algo- 
rithms to deal with the above optimisation problem. 
Later, Xiao and Dougherty [94] provided a construc- 
tive algorithm for structural intervention and applied it 
to a WNT5A network. The proposed algorithm focuses 
on the impact one-bit predictor function perturbations 
have on state transitions and attractors. Their approach, 
however, does not directly characterise the steady- 
state distribution changes that result from (structural) 
perturbations of a given probability. In order to solve this 
problem, Qian and Dougherty [95] derived a formal char- 
acterisation of optimal structural intervention, based on 
the general perturbation theory in finite Markov chains. 
Specifically, they gave an analytical solution for comput- 
ing the perturbed steady-state distribution by looking at 
function perturbations. Their work mainly focused on 
one-bit function (or rank-1 matrix) perturbations, imply- 
ing that for more general perturbations, one needs to 
consider an iterative approach. The associated complexity 
of such an approach is of 0(2 3w ), where n is the num- 
ber of genes in the network. Their results have been 
applied to a WNT5A network and a mammalian cell 
cycle related network, respectively. More recently, Qian 
et al. [96] extended their previous result in [95] to a 
more efficient solution that uses the Sherman-Morrison- 
Woodury (SMW) formula [97] to deal with rank-/c matrix 
perturbations. Thus, they managed to reduce the com- 
putational complexity of the approach from 0(2 3n ) to 
0(/c 3 ), where k 2 n (k is much smaller than 2 n ). 
The application of the derived structural intervention 
method to a mutated mammalian cell cycle network 
shows that the intervention strategy can identify the 
main targets to stop uncontrolled cell growth in the 
network. 

Qian and Dougherty [98] also looked at how long-run 
sensitivity analysis can be used in PBNs, in terms of 
difference between steady-state distributions before and 
after perturbation, and with respect to different elements 
of the network, e.g., probabilistic parameters, regulatory 
functions, etc. 

External control 

While structural intervention focuses on a permanent 
change in the network dynamics, external control relies on 
Markov decision processes theory for driving a network 
out of an undesired state, i.e. as to reach a more desirable 
one [8,18]. 

The first approach to deal with PBNs was proposed by 
Shmulevich et al. [18]. They studied the impact of ran- 
dom gene perturbations g on the long-run behavior of a 
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network. The main idea of Shmulevich et al [18] is to 
construct a formulation of the state-transition probability 
that relies on the probability of a gene perturbation and on 
Boolean functions for finding bounds for the steady-state 
probability. Their particularly interesting finding is that 
these states (which in terms of mean first-passage times 
(MFPT) are easy to reach from other states) are more sta- 
ble with respect to random gene perturbations. In gene 
regulatory networks, it is important to identify what genes 
are more likely to lead the network into a desirable state 
when perturbed. MFPT naturally captures this idea - a 
few other methods developed by Shmulevich et al. [18] 
work, for example, by maximising the probability to enter 
some particular state in some fixed maximum amount 
of time, or by minimising the time needed to reach 
that state. 

Gene perturbation works by single flips of a genes 
state, providing a natural platform for external interven- 
tion control via auxiliary input variables. It makes sense 
from a biological perspective, for example, to model aux- 
iliary treatments in cancer such as radiation. The value 
of these variables can be thus chosen such as to make 
the probabilistic distribution vector of the PBN evolve in 
some desired manner. 

More formally, given a PBN with n genes and k 
control inputs, ui, U2, . . . , u^, the vector u(t) = 
(ui(t)> U2(t), . . . , Ufr{t)) is used to denote the values of all 
control inputs at a given time step t. Let P denote the tran- 
sition probability matrix of the PBN, evolving according 
to w(t+ 1) = w(t) -P(u(t)). It is obvious to see that, at each 
time step t, P depends not only on the initial probability 
distribution vector, but also on the values of the control 
inputs. External control is essentially about making the 
network evolve in some desired manner by choosing, at 
each time step, input control values. The sequence of con- 
trol inputs, referred to as a control policy or strategy, can 
be associated to a cost function which has to be minimised 
over the entire class of allowed policies. Such functions 
capture the cost and benefit of using interventions, and 
are normally application dependent. For the sake of sim- 
plicity, we use Jco(z(0)) to denote the cost with respect to 
a control policy co and an initial state z(0). Then, an opti- 
mal PBN control problem can be defined as a search for a 
control policy co that minimises the cost /^(z(0)). External 
control in PBNs can be classified into the following two 
groups. 

Finite-horizon external control 

The finite-horizon external control problem is about mod- 
ifying over a transient period of time the network dynam- 
ics of some given PBN, without changing its steady-state 
distribution. In other words, external control is only 
applied over a finite number of M time steps, using 
policies of the form co = (/xn, /xi, . . . , /xm-i). The first 



optimal finite control formulation in PBNs, and a solu- 
tion based on Dynamic Programming [99], were given by 
Datta et al. [100]. Working assumptions implied known 
transition probabilities and horizon length, later removed 
in [101] by making use of measurements, thought to be 
related to the underlying Markov chain states of the PBN. 
Pal et al. [17] extended the results of Datta et al. [100,101] 
to context-sensitive PBNs with perturbation. The results 
have been used to devise a control strategy that 
reduces the WNT5A genes action in affecting biological 
regulation. 

Optimal finite-horizon dynamic programming based 
control, assuming a fixed number of time steps M and 
a fixed number of controls /c, has a computational com- 
plexity of 0(2 2 ), where n is the number of genes in the 
network. Namely, the problem is limited by the size of 
the network as one needs to compute the transition prob- 
ability matrix. In particular, Akutsu et al. [102] proved 
that the problem is NP-hard. h Chen and Ching [103] used 
dynamic programming in conjunction with state reduc- 
tion techniques [104,105] to find an optimal control policy 
for large PBNs. They managed to reduce the computation 
complexity to 0(| R |), where | R | is the number of states 
after state reduction. 

Kobayashi and Hiraishi [106] proposed an integer pro- 
gramming based approach that avoids computing the 
probability matrix in optimal finite-horizon control. Later, 
they extended their work to context-sensitive PBNs 
[107,108], focusing on the lower and upper bounds of the 
cost function. Furthermore, Kobayashi and Hiraishi [109] 
proposed a polynomial optimisation approach where a 
PBN is first transformed into a polynomial system, sub- 
sequently allowing to reduce the optimal control to a 
polynomial optimisation problem. In the above papers, 
only small examples are used to illustrate the proposed 
approaches. 

Ching et al. [110] looked at hard constraints for an upper 
bound on the number of controls, and proposed a novel 
approach that requires minimising the distance between 
terminal and desirable states. They also gave a method to 
reduce the computational cost of the problem by using 
an approximation technique [12]. Cong et al. [Ill] made 
one step further by considering the case of multiple hard 
constraints, i.e., the maximum numbers of times each 
control method can be applied, developing an algorithm 
capable of finding all optimal control policies. A heuris- 
tic approach was developed by the same authors in order 
to deal with large size networks [111]. A different and 
more efficient algorithm, using integer linear program- 
ming with hard constraints, was presented later by Chen 
et al. [112]. The WNT5A network is a typical example 
used in [111,112]. 

Instead of minimising the cost, Liu et al. [113] investi- 
gated the problem of how control can be used to reach 
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desirable network states, with maximal probability and 
within a certain time. Later, Liu [19] imposed another 
new criterion for the optimal design of PBN control 
policies, namely the expected average time required to 
transform undesired states into desirable ones. In both 
papers, the optimal control problem can be solved by 
minimising the MFPT of discrete-time Markov decision 
processes. 

The controllability problem of PBNs was studied by 
Li and Sun [114]. A semi-tensor product of matrices, as 
described in their work, allows to convert a probabilistic 
Boolean control network into a discrete time system. They 
provided some conditions for the controllability of PBNs 
via either open or closed loop control. 

Infinite-horizon external control 

Infinite-horizon external control implies working with 
external auxiliary variables, over an infinite period of time, 
the steady-state distribution being also changed. Policies 
in this case have the form of co = (/xn, /xi, . . .). 

In the finite-horizon case, the optimal control policy 
is calculated by (essentially) using a backward dynamic 
programming algorithm, ending once the initial state is 
reached. However, this approach cannot be applied to 
infinite-horizon control directly due to the non-existence 
of a termination state in the finite-horizon case, poten- 
tially leading to an infinite total cost. Pal et al. [115] 
extended the earlier finite-horizon results to the infinite- 
horizon case for context-sensitive PBNs. They solved 
the above two problems by using the theory of average 
expected costs and expected discounted cost criteria in 
Markov decision processes. For applications, they consid- 
ered a gene network containing the genes WNT5A, pirin, 
S100P, RET1, MARTI, HADHB, and STC2. 

A robust control policy can be found in Pal et al. 
[116], devised via a minimisation of the worst-case 
cost over the uncertainty set, with uncertainty defined 
with respect to the entries of the transition probability 
matrix. 

Due to the computational complexity of 0(2 2 ), sev- 
eral greedy algorithms have been proposed in the lit- 
erature. Vahedi et al. [117] developed a greedy control 
policy that uses MFPT. Their main idea is to reduce 
the risk of entering undesirable states by increasing (or 
decreasing) the time needed to enter such a state (or, 
respectively a desirable state). Performance of the MFPT- 
based algorithm was studied on a few synthetic PBNs 
and a PBN obtained from a melanoma gene-expression 
dataset, where the abundance of messenger RNA for 
the gene WNT5A was found to be highly discriminating 
between cells with properties associated with high or low 
metastatic competence. Later, three different greedy con- 
trol policies were proposed by Qian et al. [118], using the 
steady-state probability mass. The first one explores the 



structural information of a basin of attractors in order 
to reduce the steady-state probability mass for undesir- 
able states, while the remaining two policies regard the 
shift in the steady-state probability mass of undesirable 
states as a criterion when applying control. The identi- 
fied three policies, together with the one based on MFPT 
[117], were evaluated on a large number (around 1000) of 
randomly generated networks and a mammalian cell cycle 
network [119]. 

Some types of cancer therapies like chemotherapy, are 
given in cycles with each treatment being followed by a 
recovery period. Vahedi et al. [120] showed how an opti- 
mal cyclic control policy can be devised for PBNs. Yousefi 
et al. [121] extended the results in [120] to obtain opti- 
mal control policies for the class of cyclic therapeutic 
methods where interventions have a fixed-length dura- 
tion of effectiveness. Both of the two approaches [120,121] 
were applied to derive optimal cyclic policies to control 
the behavior of regulatory models of the mammalian cell 
cycle network [119]. While the goal of control policies 
is to reduce the steady-state probability mass of unde- 
sirable states, in practice it is also important to limit 
collateral damage, to consider when designing control 
policies. Based on this observation, Qian and Dougherty 
[122] developed two new phenotypically-constrained con- 
trol policies by investigating their effects on the long-run 
behaviour of the network. The newly proposed policies 
were examined on a reduced network of 10 nodes. The 
network was obtained from gene expression data collected 
for the study of metastatic melanoma (e.g, see [91]). 

Relationship between PBNs and other probabilistic 
graphical models 

Probabilistic graphical models, commonly applied in com- 
putational biology for network reconstruction, provide 
the means for representing complex joint distributions. 
Examples include PBNs, Bayesian networks and their vari- 
ants, e.g., dynamic and hierarchical Bayesian networks, 
hidden Markov models, factor graphs, Markov random 
fields, conditional random fields, Markov logic networks, 
etc. In this section we discuss the relationship between 
the two of them which are usually employed to deal with 
system dynamics: the PBNs and the dynamic Bayesian 
networks, the latter generalising hidden Markov models. 

A Bayesian network is essentially a graphical, com- 
pact representation of a joint probability distribution. 
The Bayesian network consists of two elements. First, a 
directed acyclic graph (DAG) where the vertices of the 
graph represent random variables and the directed edges 
or lack thereof encodes the so-called Markovian assump- 
tion, which states that each variable is independent of its 
non-descendants, given its parents [8,123]. Second, a set 
of local conditional probability distributions for each ver- 
tex, given its parents in the graph. By the chain rule of 
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probabilities, the joint probability distribution on the ran- 
dom variables in the graph can be decomposed into a 
product of the local conditional probabilities, i.e., if there 
are n random variables X/, i = 1, 2, . . . , n and Pa(X/) 
denotes the parents of X{ in the graph, then the joint 
probability distribution factors as 

n 

Pr(X!,X 2 , . . . ,X n ) = f] Pr(X;|Pa(X;)). (8) 

i=l 

Two different Bayesian networks can encode the same 
set of independencies. Such networks are said to be equiv- 
alent Equivalent networks cannot be distinguished when 
inferring the network from measurement data. One way 
to bypass this difficulty is to perform targeted interven- 
tion experiments which can narrow the range of possible 
network architectures. 

Dynamic Bayesian networks (DBNs) are extensions of 
Bayesian networks to the temporal domain and can be 
used to model stochastic processes [70]. DBNs generalise 
hidden Markov models and linear dynamical systems by 
representing the conditional dependencies and indepen- 
dencies between variables over time. Contrary to Bayesian 
networks, DBNs can be used to model feedback rela- 
tionships, a ubiquitous element in genetic regulation. In 
comparison to PBNs, dynamic Bayesian networks support 
the assignment of quantitative state values, making this 
modelling approach more flexible to handle various types 
of data. DBNs are broadly applied to represent biologi- 
cal networks such as gene regulatory networks [124-127], 
signal transduction networks, e.g., [128-130], metabolic 
networks [131], as well as networks in physiology and 
medicine [132-136]. 

As shown in [137], PBNs and binary- valued DBNs 
whose initial and transition Bayesian networks are 
assumed to have only within and between consecutive 
slice connections, respectively, can represent the same 
joint probability distribution over their common variables. 
This is true both for independent as well as dependent 
variants of PBNs. However, there are many statistically 
equivalent PBNs that correspond to a DBN. On one hand, 
the PBN framework can be considered as redundant from 
the probabilistic point of view. On the other hand, it is 
richer from the functional point of view because it models 
the regulatory roles of different gene sets in more detail 
than the conditional probabilities in DBNs [137]. The 
conversion algorithms between the two modelling for- 
malism are presented in [137], both for independent and 
dependent PBNs. Also the extensions of standard PBNs 
to context-sensitive PBNp is discussed. The perturbations 
and context switching can be introduced in the DBN for- 
malism by adding additional hidden nodes to the dynamic 
Bayesian network, as shown in [137]. 



In terms of applications, it has been shown that both 
the PBN and the DBN approaches principally have good 
performance on the inference of gene regulatory networks 
from microarray data [138]. In addition, the connection 
between PBNs and DBNs makes it possible to apply the 
advanced DBNs to PBNs tools and vice versa. For example, 
an abundant collection of learning theory and algorithms 
for DBNs already exists and methods for the analysis 
of temporal behaviour of DBNs are already established. 
These techniques can be tailored to be applied directly in 
the context of PBNs. Conversely, the tool for controlling 
the steady-state behaviour of the networks, tools for net- 
work projection, node adjunction, resolution reduction as 
well as efficient learning schemes can be applied to DBNs. 

As presented in [139], PBNs and dynamic Bayesian 
networks can be viewed as consisting of a probabilis- 
tic (Markov chain) and of a (Boolean) logic component. 
In the case of a dynamic Bayesian network, the proba- 
bilistic component is defined by a conditional probability 
chain rule and a Markov chain while the logic component 
is given by propositional logic with structural require- 
ments. As shown in [139], Bayesian networks, with their 
hierarchical and dynamic variants, as well as probabilis- 
tic Boolean networks, are all generalised by Markov logic 
networks. The same separation of components applies. 
For a Markov logic network, the probabilistic compo- 
nent is a Markov random field and the logic compo- 
nent is the first order logic. We refer to [139] for more 
details on this framework, its applications in biology and 
medicine as well as the relationship with Bayesian net- 
works. 

PBN applications in biological and biomedical 
studies 

PBN models for the representation of biological networks 

Even though a significant part of the research on PBNs is 
theoretical, a large number of applied studies on the use 
of PBNs for various biological systems can be found in 
the literature. This is particularly the case with inference 
of models for molecular and physiological networks (from 
prior knowledge or data), with subsequent model analysis 
that leads to novel knowledge in biology and medicine. 

PBNs as models of gene regulatory networks 

PBNs were originally developed as models for Gene Reg- 
ulatory Networks (GRNs) [3,8]. As stated in [32], PBNs 

1) incorporate rule-based dependencies between genes; 

2) allow the systematic study of global network dynam- 
ics; 3) are able to cope with uncertainty, both in the 
data and model selection; and 4) permit the quantifica- 
tion of the relative influence and sensitivity of genes in 
their interactions with other genes. In the PBN modelling 
framework, gene expression is quantised to two levels: ON 
and OFF. 
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The dynamical behaviour of PBNs can be used to model 
many biologically meaningful phenomena, such as cellular 
state dynamics possessing switch-like behaviour, hystere- 
sis, stability, and etc. [32,140]. Often, the attractor cycles 
are interpreted as functional states on physiological time 
scales or as cellular phenotypes on developmental time- 
scales [7,8]. This interpretation is fairly reasonable as most 
cell types are characterised by stable recurrent patterns of 
gene expression [31]. 

In the past years, there were several studies which suc- 
cessfully applied PBNs for the construction of GRNs from 
high-throughput gene expression microarray experiment 
data. In 2006, Yu et al. inferred a GRN of the inter- 
feron pathway in macrophages using time-course gene 
expression data [22]. The optimal network was identi- 
fied applying the CoD approach. It was shown that the 
respective selection probabilities are varying for different 
biological conditions, e.g., after interferon treatment or 
after viral infection on macrophage, while the structure of 
the constituent network, i.e., predictor functions, remains 
stable. With a similar approach, Nguyen et al. inferred a 
GRN of hepatocellular carcinoma from microarray data 
and compared it to a network derived from control non- 
cancerous samples [141]. They indicated that certain 
genes in tumour samples show activity in steady-state 
periods while there is no activity for these genes in the 
control (non-cancerous) samples. This allowed to distin- 
guish different gene regulatory processes being realized 
with the same set of genes. 

Hashimoto et al. modelled the cell cycle of budding 
yeast by using context-sensitive PBNs [23]. They showed 
that the switching behaviour from stationary Gl phase 
to excited Gl phase in the PBN model is more frequent, 
when compared to the stochastic model of Zhang et al. 
[142]. Recently, Todd et al. identified the ergodic sets of 
states in PBNs that correspond to each phase of the bud- 
ding yeast cell cycle, which in turn correspond to the 
cellular phenotypes [44]. The analysis of the dynamical 
behaviour gave additional insights on yeast cell cycle regu- 
lation, e.g., the yeast cell cycle network showed robustness 
both to external variable environments and to certain per- 
turbations such as nitrogen deprivation, where yeast cells 
proceeded through one round of division and arrest at Gl 
phase without appreciable growth. 

In 2011, Flottmann et al. modelled the regulatory pro- 
cesses that govern the production of induced Pluripotent 
Stem (iPS) cells by considering the interplay between gene 
expression, chromatin modification, and DNA methyla- 
tion [24] . As there is no clear guideline on how to assign 
Boolean functions to represent the interactions of each 
gene, their PBN model was designed to work by repre- 
senting uncertainty via two assignments. First, a number 
of possible functions were assigned to the corresponding 
nodes with different probabilities. Second, the influences 



of certain nodes were split into separated Boolean func- 
tions with varied selection probabilities. A flexibility was 
thus allowed for choosing Boolean functions that fit the 
experimental data. With their PBN model, an extensive 
analysis was performed, allowing to demonstrate epige- 
netic landscape changes from differentiated cells to iPS 
cells as a function of time step. In addition, by looking at 
model variants of the core iPS regulation, it was shown 
that an increased chromatin modification rate could 
improve reprogramming efficiency while faster changes in 
DNA methylation could provide an enhanced rate though 
at the price of trading-off efficiency. 

PBN within signal transduction network and metabolic 
network modelling 

To date, there is no study which specifically applied 
PBN as a stand-alone framework for modelling sig- 
nal transduction or metabolic networks. Nevertheless, 
PBN was combined with other algorithms or modelling 
frameworks. Fertig et al. presented GESSA, Graphically 
Extended Stochastic Simulation Algorithm, a mechanis- 
tic hybrid model which integrates the network model 
of cell signalling with pooled PBN to a differential 
equation-based model of transcription and translation 
computed by a stochastic simulation algorithm [25]. 
The cell signalling PBN model is generated by simu- 
lating individual protein copies with the correspond- 
ing state transitions updated according to the rules in 
the PBN. The sum of the resulting molecular states 
across copies, i.e., of each individual species, is com- 
pared to the initial state, the difference being afterwards 
returned and the cellular state being updated. GESSA 
was applied to the study of the cell fate decision of val- 
val precursor cells in C. elegans, where model predic- 
tions matched the experimental results even for mini- 
mal parameterisations of the PBN. It was thus shown 
that PBN could be an essential component when flexi- 
bility is needed in multi-level data integration and model 
construction. 

In metabolic modelling, Chandrasekaran et al. pre- 
sented an automated algorithm for the Probabilistic Reg- 
ulation of Metabolism (PROM), allowing to reconstruct 
a probabilistic GRN integrated with a metabolic net- 
work from high-throughput data[26]. PROM makes use of 
conditional probabilities to model transcriptional regula- 
tion, similar to the CoD concept in PBN inference. This 
formalism permits the strength of transcription factor 
(TF)-gene regulation as well as gene states to be rep- 
resented in terms of probabilities. PROM was used to 
generate a genome-scale integrative transcriptomic and 
metabolomic network of Escherichia coli, where PROM 
surpassed the state-of-the-art methods such as the regu- 
latory flux balance analysis. PROM was also used to gen- 
erate an integrative model of Mycobacterium tuberculosis. 
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The results from the model analysis offered additional 
details on known regulatory mechanisms and also helped 
to uncover the function of less studied genes on metabolic 
regulation. 

Apart from these two studies, several other works also 
made use of a probabilistic framework for analysing sig- 
nal transduction and metabolic networks. Kaderali et al, 
for instance, developed an algorithm that reconstructs 
signalling pathways from gene knockdown data (RNAi 
data) [143]. In this work, pathway topologies are inferred 
by using Bayesian networks with probabilistic Boolean 
threshold functions. The algorithm was used to study 
the Janus Kinase and Signal Transducers and Activators 
of Transcription (JAK/STAT) pathway, correctly recon- 
structing the core topology of the pathway along with 
model variants. Similarly, Sauer et al. [144] used prob- 
abilistic equations to determine flux ratios, allowing to 
express the relative contribution of certain metabolites or 
pathways as modulators in the network. This assignment 
is more realistic than using flux absolute integer numbers, 
given that the flux of each source can relatively contribute 
to the production of certain metabolites. 

PBN applications in the context of physiology 

PBNs were also used in the recent years for studying net- 
works in physiology, with a close link to medicine. Tay 
et al. described a dengue hemorrhagic fever (DHF) infec- 
tion model which contains the interplay between dengue 
virus and different cytokines which are cross-regulated in 
T-helper 1 (Thl) and Th2 cells [9]. In their work, a sin- 
gle probabilistic Karnaugh-Map is generated, modelling 
the inducement probability of each cell as to define the 
overall influence of inducing nodes. Simulation results 
matched clinical data for both synchronous and asyn- 
chronous updating, with respect to the form and the 
average duration-based attractors, respectively. In addi- 
tion, by applying a genetic algorithm [145] to modulate 
the DHF attractor basins to dengue fever (DF) basins (a 
less severe form of DHF), Tay et al. also identified the 
tumour growth factor beta (TGF/3), interleukine-8 (IL-8) 
and IL-13, as sensitive intervention points. 

Another example in this field can be found in the study 
of Ma et al, where, based on functional Magnetic Reso- 
nance Imaging (fMRI) data, the authors developed a brain 
connectivity network model for Parkinsons disease [10]. 
A method similar to the one of Yu et al. [22] was used 
for probability inference selection, i.e., the calculation of 
CoDs. Then the CoDs were subsequently used to gener- 
ate an influence matrix representation of the brain sig- 
nal connectivity among brain components. The obtained 
results showed that a significant difference in connec- 
tivity exists for many paired brain-components com- 
paring between normal, Parkinsons disease with drug, 
and Parkinsons disease with drug withdrawal conditions, 



and this difference was expressed in terms of estimated 
range of coefficient mean activity. This particular infor- 
mation may allow to construct a new screening procedure 
for Parkinsons disease diagnose and to determine drug 
trial responsiveness based on a non-invasive, fMRI-based 
investigation in the future. 

A certain number of the previously described (applied 
research) articles on PBN have applications not only in 
molecular biology, but also in physiology or medicine. 
Only to name a few examples, being able to distinguish 
among the regulatory networks of cancer and healthy 
cells, as presented by Nguyen et al, could contribute to an 
early detection of cancerous genes in susceptible popula- 
tions [141]. A better understanding of dynamic processes 
and the control of somatic cell programming, as proposed 
by Fertig et al, may lead to a future use of iPS cells in cell 
or tissue replacement therapies [25]. Last but not least, 
the PROM algorithm, as introduced by Chandrasekaran 
et al, is capable of predicting transcription factor drug 
targets which are major hubs in the cellular network of 
pathogenic organism such as Mycobacterium tuberculosis 
[26]. A further development of drugs in this direction may 
help in the treatment of different infectious diseases. This 
new line of treatment could have a strong impact for third- 
world countries where infectious diseases still remain a 
major cause of death. 

PBN for Systems Biology and Systems Biomedicine? 

As previously discussed, the PBN framework is a topic of 
intensive and continuous theoretical research with suc- 
cessful applications in the biomedical area. To describe 
and extend a vision on future PBNs' applications, we sum- 
marise additional arguments to support why this mod- 
elling approach is suitable for future research in Systems 
Biology and Systems Biomedicine. 

Data integration 

Different types of biological and clinical investigation 
datasets, ranging from qualitative to high-throughput 
quantitative experimental data, were successfully applied 
in PBN inference and analysis. Yu et al. [22] and Nyugen et 
al. [141], for instance, inferred GRNs of macrophages and 
hepatocellular carcinoma using microarray gene expres- 
sion data. Flottmann et al. [24] built a comprehensive 
epigenetic regulatory network of iPS cells based on gene 
expression, chromatin modification and DNA methyla- 
tion data generated from multiple high- throughput exper- 
iments. Ma et al. applied voxel selection on fMRI clinical 
data to capture the activities of each brains compartment 
as the inputs for learning a functional brain connectivity 
network [10]. 

We have recently shown that the normalised activity of 
signalling proteins from quantitative western blot exper- 
iment can be compared to the steady-state probability 
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of certain molecule to be ON in instantaneously-random 
PBNs. In an ergodic model, the activities of signalling 
proteins, usually given by their phosphorylated forms nor- 
malised to the maximal signal, could be correlated with 
the steady-state probability distribution on the state space 
of the PBN model. With this regard, PBN could support 
the integration of semi-quantitative experimental data. 
Apart from quantitative western blot data, the profiles 
of signalling proteins from alternative experiments such 
as enzyme-link immunosorbent assay (ELISA) and high- 
throughput protein array data are also compatible with 
this framework (publication submitted). 

The PBN framework also allows for the description 
and analysis of large-scale models, for instance as in 
the case of a Boolean model of apoptosis of Schlatter 
et al. [146]. Therein, a PBN model was derived from 
the original literature-based BN consisting of 86 nodes 
and 125 Boolean interactions. Quantitative experimental 
data in this study were normalised to the maximal sig- 
nals across experiments and were used as input data for 
the PBN model. We analysed the strengths of canoni- 
cal pathways and crosstalk interactions between different 
signalling components among apoptotic and related sig- 
nalling pathways through the identification of selection 



probability. It was possible to obtain these via optimi- 
sation. Thereby a curated signal transduction network 
topology was derived. The resulting PBN demonstrates 
the correlation between UVB irradiation, NF/cB, caspase 
3, and apoptotic activities in a semi-quantitative manner 
which could not be demonstrated by the original BN. The 
analysis pointed at an inconsistent caspase 3 measure- 
ment, which shows no activity for high UVB irradiation 
while significant apoptosis is measured (see Figure 6, 
publication submitted). 

Furthermore, the PBN framework has a good poten- 
tial to describe cellular dynamics at multiple levels. 
Hybrid PBN-related models could be applied, as previ- 
ously described, e.g., in the studies of Fertig et al. and 
Chandrasekaran et al. [25,26]. As reviewed in detail by 
Goncalves et al. [147], bridging layers towards an inte- 
gration of signal transduction, regulation and metabolism 
into mathematical models still posts many challenges as 
each of the biological layer has their own distinct char- 
acteristics and therefore is suitable for only a subset of 
modelling approaches. To address such challenges, an 
integrative hybrid model for flux balance analysis was pro- 
posed, combining BN modelling for the gene regulatory 
part, ODE modelling for the signal transduction part and 
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flux balance analysis for the metabolic part. With this 
regard, PBN could also be integrated as part of such a 
hybrid model to describe GRNs and/or signalling net- 
works to provide more details on modelling analysis and 
interpretation comparing to traditional BNs. 

Computational tools for PBN modelling and analysis 

Several PBN modelling and analysis tools were continu- 
ously developed over the past recent years. The BN/PBN 
MATLAB-based toolbox, introduced by Lahdesmaki and 
Shmulevich in 2003 [27], deals with the simulation, analy- 
sis (network statistics, state transitions and distributions), 
visualisation and intervention analysis of both BN and 
PBN models. The toolbox was specifically designed for 
GRN inference and it makes use of CoD calculations. 
State transition probabilities and influence values (the 
indicators for interactive effect for each pair of genes) 
are subsequently calculated based on these calculated 
CoDs. Ma et al. successfully applied the BN/PBN tool- 
box to infer and analyse the brain connectivity network 
of Parkinsons disease patients, as previously described in 
Section 'PBN applications in the context of physiology! 

Hinkelmann et al. introduced ADAM (Analysis of 
Discrete Models of biological systems using computer 
Algebra) [28], a web-based tool for rapid steady-state 
identification in various discrete model types. The tool 
automatically converts discrete models into polynomial 
dynamic systems, allowing to run computer-based alge- 
bra analysis. For probabilistic networks, ADAM generates 
a graph of all possible (local rule) updates, thus being 
capable to build an enumeration of all steady states. Bool- 
net, as introduced by Mussel et al, is an R-package for 
the generation, modelling, reconstruction and analysis 
of both synchronous and asynchronous BNs or PBNs 
[29]. The toolbox features time-series (experimental data) 
based network inference, e.g., making use of Markov chain 
simulations for attractor identification with subsequent 
visualisation and robustness analysis via network pertur- 
bation or heuristic search and random walks. We have 
recently developed optPBN, a MATLAB-based toolbox 
for PBN optimisation based on the BN/PBN toolbox. 
PBNs can easily be constructed from Boolean rule-based 
models. The toolbox also provides a flexible platform 
for data integration (e.g., to integrate data from multi- 
ple experiments). Different algorithms can be used to 
address the resulting optimisation problem. Thus, based 
on normalised protein activity at steady-state data, one 
can identify a curated model structure from different can- 
didate models. Subsequent analysis on the curated PBN 
can be performed in the BN/PBN toolbox (publication 
submitted). 

We also discuss a few different algorithms and tools 
which are not specifically designed for PBN but with 
a high potential for the analysis of PBNs. PROM, for 



example, offers a mean to calculate the flux activities 
of a metabolic network in a probabilistic manner based 
on gene expression data [26]. Specifically, this gives rise 
to the applicability of the PBN framework for metabolic 
models. Recently, Terfve et al. introduced CellNOptR, a 
flexible toolkit for training protein signalling networks 
based on a multiple logic formalism [148]. CellNOptR 
offers support for optimisation with respect to multiple 
modelling frameworks, ranging from logical to ODE (logic 
rule derived) models. Extending CellNOptR towards a 
probabilistic modelling framework is also foreseen for 
future work. 

A perspective on potential applications of PBNs in a clinical 
setting 

It has been a decade since the completion of the Human 
Genome Project in 2003 that initiated the era of bio- 
logical and medical investigation in omic scales [149]. 
Due to technology advancements, the costs of genome 
sequencing and high throughput biomedical investiga- 
tions are exponentially decreased and they might become 
part of the routine medical investigations in a fore- 
seeable time frame [150]. Datasets from omics exper- 
iments usually consist of large lists of numbers that 
represent genes, transcripts, proteins, or metabolites 
depending on the method applied. In the near future 
all these methodologies might be applied together rou- 
tinely, even in time series examinations. The major 
problem with such data is their high complexity and 
the need to make them interpretable by the medical 
staff. Therefore, there is a strong demand for reason- 
able computational approaches to integrate multidimen- 
sional "big data" [151]. In addition, given the rich sets 
of information from individual patients that physicians 
will acquire, smart approaches are mandatory to trans- 
late and simplify these large-scale biomedical data. Such 
approaches should facilitate a physicians decision-making 
process to provide more accurate diagnosis and optimal 
treatment. 

For these fields we identify the PBN framework as 
a powerful tool. Recent applications of PBN modelling 
of gene regulatory and signalling networks have been 
described in the previous section. As previously sum- 
marised, PBNs allow an effective visualisation of GRN 
models [9,10], allowing to represent gene function and 
activity [152]. These efforts foster the understanding of 
gene-gene interactions, consequences of aberrant gene 
function and targeted perturbations of such networks, 
as well as finding out the least adverse effects of per- 
turbations [9,153]. PBNs allow for the integration of 
information from large data sets and for inferring log- 
ical relationships between genes/networks. This feature 
is of particular benefit as many relationships and struc- 
tural connections among genes are not known. Unknown 
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relationships between transcripts and proteins can also 
be assessed. In a therapeutic perspective PBNs could be 
used in a disease-relevant context because many, fore- 
most chronic diseases, share probably common under- 
lying mechanisms that are not elucidated so far [154]. 
Using PBNs in the study of disease-related networks could 
enable us to take genetic interactions into account and 
associations could be generated to identify comorbidities 
sharing common causative factors. Skahanenko et al. for 
instance have applied Markov logic networks, a proba- 
bilistic logic modelling approach in the same category as 
PBN, to explore gene-phenotype associations. Whereas 
traditional statistical methods are employed to identify 
the marker that associates the most with an observed phe- 
notype, Markov logic networks can be used to identify a 
subset of markers that predicts the phenotype. Within this 
method, the relationship between the genetic markers and 
phenotype(s) can be hypothesised and modelled. All mod- 
els can then be tested and their respective probability can 
be derived [139]. 

In the context of a single, yet complex disease, the 
study on brain connectivity in Parkinsons disease by Ma 
et al. [10] is a good example showing how a probabilistic 
model such as PBN could translate large-scale biomedical 
data into a potential application in clinic. fMRI princi- 
pally measures blood oxygen level-dependent signals that 
are correlated to the blood flow into different regions 
of the brain, which in turn give physicians information 
on the functional activity of specific areas of the brain 
[155]. For some neurological disorders, such as Parkin- 
sons disease, the lesions mainly affect a specific area of 
the brain such as basal ganglia, but have consequences 
on the overall integrity of brain connectivity, especially 
on the dopaminergic pathway-dependent motor and cog- 
nitive control [156]. Therefore, considering the aetiology 
and disease progression from only conventional MRI data 
which demonstrate only structural information is cer- 
tainly insufficient to yield a comprehensive understanding 
on the course of disease. Considering diseases as net- 
work perturbations [157], the PBN model from Ma et al. 
demonstrated differences of brain connectivity networks 
comparing healthy population and diseased cases with 
and without medication. Such observations could pos- 
sibly be further developed towards clinical biomarkers 
which could then be added to physicians' portfolio and in 
turn facilitate diagnostic process, treatment design, and 
follow-up strategy. 

Generally, the incorporation of tentativeness and prob- 
ability could be evolved into a valid concept in a clinical 
setting, as routine medical investigation often provides no 
conclusive data. Together with a comprehensive reduc- 
tion and translation of large-scale and complex biomed- 
ical data, the PBN framework might serve as a mean to 
develop simplistic terms like a probability score for certain 



condition, e.g., for having a disease or of being respon- 
sive to treatment. Such a probabilistic score could serve 
as a simple but powerful additional input for physicians 
in order to improve their healthcare management. As a 
whole, healthcare systems would benefit from reducing 
costs related to unnecessary diagnostic investigations and 
treatment failures. 

Conclusion 

Even though the concept of PBN for the modelling of 
biological systems is still young compared to other mod- 
elling approaches, a broad area of research activities on 
this modelling approach such as network inference and 
network control have been well-established and are con- 
tinuously developed. For a meaningful comparison of dif- 
ferent inference algorithms in the future, it is necessary 
to quantify their performance. The prospective research 
in the area of network inference is to develop a formal 
framework for validation of network inference proce- 
dures. Moreover, there is a demand for establishing the 
properties of network inference procedures under vari- 
ous conditions, e.g., model class, distance function, etc. 
The current trend in structural intervention and exter- 
nal control is to develop new methods to reduce their 
computational complexity and to define the optimal con- 
trol problems and find the corresponding optimal policies 
for specific therapies. With its flexibility for data inte- 
gration and the availability of supporting algorithms and 
computational tools, PBN is one of the most suitable 
modelling frameworks to describe and analyse complex 
biological systems from molecular to physiological levels 
with possible future application at clinical level. 

Endnotes 

a In general, yi, /2> • • • > Yn need not be independent and 
identically distributed random variables, but for the 
simplicity of presentation are assumed so. 

b A state in a Markov chain is said to be ergodic if 
returns to the state can occur at irregular times and the 
state is positive recurrent. If all states in an (irreducible) 
Markov chain are ergodic, then the chain itself is said to 
be ergodic. 

c In a generalised PBN framework a network variable 
can have any value in {0, 1, . . . , d — 1}, where d > 2. 

d In the graph-theoretical terminology the notion of an 
ergodic set of states in a Markov chain corresponds to the 
notion of a bottom strongly connected component in a 
graph. 

e In computer science, the complexity of a function or 
an algorithm is expressed or characterised using the big 
O notation, namely, how the function or algorithm 
responds to changes in its input size. 

Optimisation deals with a broad range of problems, 
relying on, for example, convex programming, optimal 



Trairatphisan et al. Cell Communication and Signaling 201 3, 1 1 :46 
http://www.biosignaling.eom/content/1 1 /I /46 



Page 22 of 25 



control, combinatorial optimisation or evolutionary 
computation paradigms; examples and additional 
information can be found by referring to [158-165] 

g A one-time gene perturbation changes the value of 
one or more genes without modifying the rules or 
probabilistic parameters of the network. 

h In computational complexity theory, NP-hard is 
a class of problems that are at least as hard as the hardest 
problems in NP (nondeterministic polynomial time). 
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